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We introduce a simple model of population dynamics which considers birth and death rates for 
every individual that depend on the number of particles in its neighborhood. The model shows 
an inhomogeneous quasistationary pattern with many different clusters of particles. We derive the 
equation for the macroscopic density of particles, perform a linear stability analysis on it, and show 
that there is a finite-wavelength instability leading to pattern formation. This is the responsible for 
the approximate periodicity with which the clusters of particles arrange in the microscopic model. 
In addition, we consider the population when immersed in a fluid medium and analyze the influence 
of advection on global properties of the model. 
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I. INTRODUCTION 

Interacting particle systems are useful models to un- 
derstand a variety of effects in fields as diverse as con- 
densed matter physics, chemical kinetics, population bi- 
ology (where they are called individual based models) or 
sociology (agent based models) . As one of the simplest 
examples one can consider an ensemble of Brownian par- 
ticles, each one dying or duplicating with given probabil- 
ities per unit of time. Several authors [1111111 have 
considered such Brownian Bug (BB) model in the context 
of population dynamics (in particular to address plank- 
ton distributions and patchiness), in the case in which 
the probabilities of death and reproduction are equal. 
Aggregation of the particles in a decreasing number of 
clusters occurs. This clustering is somehow surprising 
since a standard mean-field or rate-equation description 
gives for the particle density p the equation 
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where D is the diffusion coefficient, and Xq and (3q are 
the birth and death rates, respectively. Obviously, when 
•^0 = Eq. is simply the diffusion equation which 
cannot lead to spatial inhomogeneities. 

This result was known since some time ago for this 
and related modelsH US H S 0, II , and points out 
the relevance of the fluctuations present in the discrete 
stochastic particle model, neglected in a nai'f mean-field 
macroscopic description, and that lead to reproductive 
pair correlations: the mean rates of death and birth are 
equal, but if locally there is an excess of reproduction 
events, and diffusion is not fast enough, a cluster of par- 
ticles will develop, whereas no birth will occur in empty 
zones and particles will simply disappear from regions 
with excess of death. 

The authors of ^] go beyond that result, and show 
that the clustering persists even in the presence of rather 
strong stirring, as it would occur if the bugs live in a 



turbulent fluid such as the Ocean. The clusters now be- 
come elongated filaments, but there is still strong spatial 
inhomogeneity arising from the microscopic particle fluc- 
tuations and reproductive correlations. 

The simple model just described misses some impor- 
tant features present in real biological populations. The 
most obvious is the absence of any interaction between 
the bugs. Among other consequences, the global dynam- 
ics of the system, i.e. the time evolution of the total num- 
ber of particles, is completely independent of its spatial 
distribution. Thus stirring the system alters the spatial 
pattern of the bugs, but neither their individual lifetimes, 
nor the time history of the particle number, nor its sta- 
tistical properties. 

In this Paper we introduce interacting particle mod- 
els by modifying the birth and death rates of the BB 
model. They will take into account the number of neigh- 
bors within a given distance of each bug. There is now 
a strong interplay between the bug dynamics and the 
ambient flow and, in addition, new effects arising from 
the spatial range of interaction occur and modify the 
reproductive-correlations clustering effect. In particular, 
an inhomogeneous steady structure with many different 
clusters of particles coming from different families (i.e. 
they are born from a different parent), and arranged in 
a periodic pattern, may occur. The number of particles 
in any of these clusters is similar, resembling the spread- 
ing of individuals in small groups over a geographical 
area. This pattern formation phenomenon occurs via a 
finite wavelength instability that can be characterized in 
a deterministic description, being ffuctuations only of sec- 
ondary importance. We analyze the phenomenon with a 
continuous-field Langevin description obtained from the 
particle model by Fock space techniques. 

The Paper is organized as follows. In the next Section 
we introduce the discrete models. In Sect. IIIII we study 
numerically some of their properties. In II VI we study the 
pattern formation process and perform a stability analy- 
sis within the continuum-field description of our model. 
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In Section [V] we study the influence of a fluid flow on the 
particle system, and in the last Section we summarize 
our Conclusions. 

Derivation of the continuum-field representation for 
our particle model is interesting by itself. In this Paper 
we will use (Section IIV|I just the deterministic part of 
this representation, but we provide a detailed derivation 
of the full Langevin equation, via Fock space techniques, 
in the Appendix, so that the additional results presented 
there can be used for future reference. Representation of 
advection processes within the Fock space formalism is 
also considered there. 



II. MODELS 

In this section we introduce the discrete models subject 
of study of the paper. We begin defining the original BB 
model and then our extensions. 



A. BB model 

The microscopic rules are simply enumerated 0- 
Let N{t) the number of bugs in the system (a two- 
dimensional periodic box of size LxL; in all our computer 
simulations we will take L = 1): 

1. There is an initial population of iV(i = 0) = Nq 
bugs or particles, randomly located. 

2. One particle is selected at random and it dies with 
probability p, reproduces with probability g, or re- 
mains unchanged with probability r (p + q + r = 1). 
In the case of reproduction, the newborn particle 
is located at the same place as the parent particle. 
The process is repeated a number N(t) of times 

lio]. 

3. Each particle moves in random direction a distance 
drawn from a Gaussian distribution of standard de- 
viation a (this models Brownian motion). 

4. When advection is considered, the particles are 
transported by an external fiow to be described 
later. 

5. Time is incremented an amount r, and the algo- 
rithm repeats. 

Symbolically, in chemical reaction notation: 

A 0, (2) 
A ^ A + A (3) 

where A represents individual particles, each one dying 
at a rate /3o = p/t (death rate per particle and unit of 
time), or reproducing at a rate Aq = q/r. The Brownian 
motion step leads to diffusion with a diffusion coefficient 



D = (t^/2t. In the following we measure time in units 
of T, SO that r = 1, /3o = p, Ao = q, and a — \J2D. We 
take r = so that Aq -I- /3o = 1 , and define the important 
parameter /i = Aq — /3o, the difference between birth and 
death rates. 



B. Neighbourhood-dependent (ND) model 

The new model is analogous to the one before, except 
that in step 2 the reproduction and death rates of a given 
particle labelled j, A(j) and /3(j), and are not constant 
but depend on the number of particles surrounding the 
particle j. Explicitly we take (with r = 1): 

m = max (^0, Ao - ^N^^ , (4) 

and 

/3(j)=max(^0,/?o-£iV^) . (5) 

where denotes the total number of particles which are 
at a distance smaller than i? from particle j (excluding 
the particle j itself). R is thus a range of interaction, 
is a saturation parameter, and a controls the asymme- 
try between its influence on death and on reproduction. 
The BB model is recovered when i? ^ 0. The maxi- 
mum condition is imposed to insure that A(j) and /3(j) 
are positive definite, as it should be given that they are 
probabilities. When Ns is positive, the model penahzes 
reproduction when particles are crowdedly surrounded. 
Lonely particles reproduce with higher probability. This 
kind of interaction would be appropriate to model in- 
dividuals that compete for resources (e.g. food), in a 
neighborhood of its actual position. When is positive, 
death rate decreases in crowded environments, modelling 
a kind of mutual protection. The opposite behavior oc- 
curs when these parameters are negative, a situation that 
will not be considered in the present Paper, though most 
of the results presented can be extended easily to this 
case. This model of interactin g p articles is related to 
many others (see reviews in 0, ) in which some lim- 
itation in the growth of the population at a single site 
is imposed via a fermionic restriction (explicitly stated 
|12l | or implicitly imposed on computer simulations by 
forbidding double occupation of lattice sites) or via the 
inclusion of the coagulation process A + A^A inverse to 
© Ji2.i 2M ■ Our model (and the BB model) shares with 
them several qualitative features, being the most impor- 
tant the fact that the empty state is an absorbing state: 
if at some moment all the population becomes extinct, no 
recovery is possible within the rules of the model. This 
leads to an absorbing phase transition from an active or 
surviving phase to an absorbing dead phase when some 
effective reproduction rate is reduced. The peculiarity 
in our model is that interactions are not purely local. 



but extend to a finite distance R. This should be ir- 
relevant for the critical behavior close to the absorbing 
phase transition, since only asymptotically large scales 
are important there, and we expect this transition in our 
model to be in the standard Directed Percolation (DP) 
or Reggeon field-theory universality class to which many 
of these interacting models belong We will see, 

however, that the behavior in the active phase is greatly 
influenced by the existence of a finite interaction range 
R. In consequence we will not analyze in great detail the 
absorbing phase transition, but concentrate in the active 
phase(s), where more novel behavior occurs. 



III. NUMERICAL STUDY OF THE DISCRETE 
MODELS 

This section is devoted to present some numerical re- 
sults that stress the differences between the BB and the 
ND model. Here we consider the system with no external 
driving flow, whose analysis is left to section Ivl 

The BB model has been studied in detail 

nasi If 

/i = Aq — /3o > 0, the total population generally explodes 
exponentially, with a time scale given by although 
there is a finite probability for extinction that depends 
on the initial population and decreases with increasing 
fi. If fj, < the final state is, with probability 1, the 
total extinction of the particles, occurring again at a ex- 
ponential rate characterized on average by fJ,~^, but with 
diverging relative fluctuations. A critical situation occurs 
when Aq = Po- In this case, the particles arising from the 
same ancestor form clusters, with the number of clusters 
decreasing in time and the number of particles in the 
surviving clusters growing, so that {N{t)) = A^o Vt, with 
the average taken over different realizations. But fluc- 
tuations in N{t) are huge (its variance diverges linearly 
in time), with some runs leading to fast extinctions, and 
others with clusters surviving for long time. In a finite 
system all clusters finally disappear, but the typical life- 
time diverges linearly with Nq, and the average life time 
is infinity f3|. Fig. ^shows the distribution of particles at 
two different stages of the evolution, one in which a large 
single cluster, coming from a single ancestor, is present 
(right panel) corresponding to a long-time evolution, and 
another with still many clusters from different ancestors 
(left panel) for an earlier stage of the temporal evolution. 

Figure [21i) shows the time evolution of the total num- 
ber of bugs N{t) in the critical case, fJ, — 0, for a partic- 
ular realization, displaying the critical fluctuations, and 
examples of cases with nonvanishing fi. One can observe 
the fast decay (growth) of N{t) for /i negative (positive). 

The behavior of the ND model is rather different. Just 
for simplicity we consider here (and in the rest of the Pa- 
per) the value a = so that only reproduction depends 
on the neighborhood. Figure I^Jd) shows the time evolu- 
tion of the population. For /i smaller than a critical value 
/ic > (which turns out to be Hc ~ 0.4 for the parame- 
ter values used in the Figure) , we find always extinction. 
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FIG. 1: Spatial configurations for the BB model at two dif- 
ferent times. Left: configuration after 100 steps, with a large 
number of surviving clusters. Right correspond to a single 
cluster remaining after 3000 steps. The value of the parame- 
ters are Aq = /3o = 0.5, D = 10~^, and the initial population 
is of A'^o = 1500 bugs randomly distributed. 



10000 



8000 



6000 



4000 



2000 



2000 r 




1000 2000 3000 1000 2000 3000 



FIG. 2: a) Total number of particles, A/'(t), vs time for the BB 
model and three different values of the control parameter 
From top to bottom: /i = 5 10~*, /j, = 0, and /i — —5 10"*; 
here D = 10~^. b) Idem for the ND model and four different 
values of /i: From top to bottom, jj, = 0.7, ^ — 0.5, n — 0.4, 
and ^ = 0.3. Two are above critical {fic ~ 0.4) and two below 
it. The other parameter values are R = 0.1, Ns — 50, and 
D = 10"^ 



whereas typical realizations reach a finite average popu- 
lation at long times for fi > fic- We plot in Fig. |31the 
total average number of particles N{t) at long-time vs 
/i, and different values of the parameters. The scaling 
used to present N{t) is suggested from an analytical ex- 
pression discussed in next Section. As discussed later, 
it provides good data collapse in the left plot Fig. |21 for 
a diffusion coefficient of D = 10^^, but it is grossly in- 
adequate for data in the right panel, corresponding to 
smaller diffusivity, D = 10"^. 

The nature of the spatial distribution in the active 
phase depends on the values of the parameters. For large 
enough D, the spatial distribution of particles is homo- 
geneous on average, whereas clear clustering occurs for 
small D. As in the BB model, the clusters are coming 
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FIG. 3: Long-time average number of particles N{t) vs fi. 
Left panel corresponds to D = 10~* and right to D = 10~^, 
and the other parameters as indicated. The average is taken 
from the instantaneous particle numbers at times between 
1000 and 10000 steps; the error bar indicates the standard 
deviation of the instantaneous fluctuations around this mean 
value. 



from different families. But here they are not transient 
and the most striking feature is that they organize in a 
periodic pattern. The periodicity of the pattern is of the 
order of R, the interaction range. In addition to decreas- 
ing D, this transition to a periodic organization occurs 
by increasing R and, for small enough D, by increasing /i. 
Fig. 21 shows examples of the different spatial patterns. 

We see that the most notable effect of the introduc- 
tion of interactions with a characteristic spatial scale has 
been the segregation of bugs in a periodic array of clusters 
(bottom- right plot in Fig. 0)). This seem to be a rather 
natural way to make compatible the high local growth 
at relatively large value of /i, with the reduction of this 
growth that a too crowded neighborhood would imply: 
the empty space between the clusters acts as a buffer 
zone keeping the competition for resources less limiting 
than in a homogeneous distribution. We expect that this 
mechanism will appear in Nature when there is a scatter- 
ing of the total population in small groups over a large 
spatial area. One can think, for instance, of the spread- 
ing of groups of predators, or even of primitive human 
societies that are aggregated in small tribes. 

We characterize the patterns in terms of the structure 
factor S{k). It is defined as 
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where Xj(i) = {xj,yj) is the position vector of the par- 
ticle j at time t, q — (q,j;,qy) is a two-dimensional 
wavevector, the sum is over all particles, and the average 
is a spherical average over all wavevectors of modulus 
|q| = K, and a further temporal average in the long- 
time state is added to improve statistics. Maxima in this 
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FIG. 4: Long-time spatial structures for the ND model. Left 
column corresponds to two patterns with the same value of 
D = 10~*, and two different values of /i = 0.5 (up) and 
fj, = 0.9 (bottom). Right column corresponds to fixed /i = 0.7, 
and D = 10"" (upper), and D = (bottom). In all the 

plots, Ns = 50 and R = 0.1. 



function identify relevant periodicities in the interparticle 
distribution. 



In Fig. |31 we show the structure factor in the steady 
state of the model for different values of the parame- 
ters. The emergence of the periodic patterns is indicated 
by the peak in the structure factor. As in Fig. 0] we 
show here two different scenarios: The upper panels cor- 
respond to the structure factor (left) of different patterns 
with pL fixed and changing D. One can observe that by 
decreasing D the value of the peak increases (upper-right 
panel) , indicating that clustering with a strong periodic- 
ity develops. Bottom panels are for D = 10~^ and differ- 
ent values of /j,. Here the pattern is rather developed at 
all values of /i above the absorbing transition at /Jc ~ 0.4, 
with only mild variations of the peak height (right panel) 
with fi. By analyzing non-spherically-averaged versions 
of (O, we confirm that the periodic pattern has hexago- 
nal symmetry at onset, although this symmetry changes 
at high values of fi. 



We next try to explain quantitatively the observed pat- 
terns in terms of an analytical description. 
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FIG. 5: Structure factors (left panels) and corresponding 
height of the main peak (right) for different patterns in the 
ND model. Upper panels are for fixed value of = 0.7 and 
the values of D shown in the legend box. Bottom panels are 
for D = lO"'^ and different /i's as shown in the legend box. In 
all the plots, Ns = 50 and R = 0.1. Note the different scales 
in all the plots. 



IV. FIELD THEORY DESCRIPTION, 
STABILITY ANALYSIS, AND SPATIAL 
PATTERNS 

Standard theories and tools for pattern formation stud- 
ies [T^ address continuous field models, and are not 
particulary well suited to analyze pattern formation in 
particle systems. Fortunately, there are well established 
techniques (known under the name of Doi-Peliti theory 
or Fock space techniques H, 0, 0) that allow a de- 
scription of interacting particle systems in terms of field- 
theoretic Langevin equations. These techniques turn out 
to be equivalent to the Poisson representation 0,1131 • In 
the simplest cases, the interacting particles are instanta- 
neously Poisson distributed in each small space region, 
and the field description gives the space and time vary- 
ing average value (/)(x, t) of the local Poisson distribution 
for the particle density p. 

In general, however, the continuous field is complex 
and this simple interpretation does not hold, but still in 
this case all the moments of the particle density p can 
be obtained from the moments of the fluctuating field 
(j). For example the first moments of both quantities are 
equal < p{x,t) >=< 0(x, t) > . In the Appendix, we 
derive in detail the Langevin field description for the ND 
model (including hydrodynamic flow, arbitrary values of 
a, etc.). Two approximations are needed to arrive to 
the final form (|46|) (see the Appendix for details). As a 
first attack to the problem of pattern formation in our 
particle model, we analyze in this Section just the deter- 
ministic part of the field equation, i.e., the noise term 
will be neglected. We will see that this will be sufficient 
to understand the main qualitative features of the pat- 
tern forming instability. The expected influence of the 
noise would be to affect system properties in the vicinity 



of transitions and instabilities, and to shift the position 
of the transition lines |15j. In addition, since our sys- 
tem is translational invariant and two-dimensional, it is 
very likely that the sharp bifurcation to patterns that we 
find in the deterministic analysis will be blurred by noise 
into a non-sharp crossover even in the thermodynamic 
limit. Nevertheless, we find good qualitative agreement 
between most of the observed properties of the discrete 
model presented in the preceding section and the deter- 
ministic predictions obtained here. 

Thus we analyze the deterministic version of Eq. 1)46(1 
(no flow and a = 0): 

dtcj){^,t)=DVm^,t) + 

(Ao - /3o)0(x, t) - ^(/.(x, t) /|,_r|</? '^(i-' ^) (7) 

At this level of mean-field- like approximation (no noise), 
the field 0(x, t) can be interpreted as the density field 
p(x, t). Stationary homogenous solutions of this equation 
are the absorbing phase 0(x, t) = 0, and the active or 
survival phase ^(x, t) = (f>s — pNa/nR'^ (remember that 
/i = Ao — /3o)- For p < the only stable solution is 
the absorbing one; the transition to the survival state is 
approached at p = 0, and this state is stable for a range 
of positive values of p. At the deterministic level the 
transition is transcritical, but we expect it to be modified 
by noise into a transition of DP type [Til ITHI | , occurring 
at values of p larger than zero, as observed. 

We make a stability analysis of the (j>s solution by 
considering small harmonic perturbations around it, 
0(x, t) — + S(f>{x,t), with 6(j){x,t) oc exp(At -I- ik • x). 
After simple calculations one arrives to the following dis- 
persion relation 

2p 



\{K) 



-DK^ 



KR 



MKR), 



(8) 



where K is the modulus of k, and Ji is the first-order 
Bessel function. It is clear that the relevant parameters 
in the problem are p and D / R? (in fact the precise adi- 
mensional combinations are pr and Dt/R^, but remem- 
ber that we are measuring times in units of r, so that 
T = 1). The eigenvalue \{K) (which is in fact a function 
of KR, p, and D/R^) is real and can be positive for some 
values of the parameters. This is shown in Fig. where 
we plot A against K for different values of p around pp 
as given below in Eq. lO, with fixed D/i?^. 
The equations 



dXjK) 
dK 
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A(X„ 



(9) 
(10) 



identify the values of the parameters at which the max- 
imum of the curve A(i^r), at Km, becomes positive. 
This gives a line of instability pp ~ pp{D/R^) in the 
parameter plane. It is straightforward to obtain that 
Pp = —DRK^/{2Ji{K„iR)) and the equation for Km 
reads 



KmR 



2MKmR) 



(MKmR) - J2{KmR)) -3 = 0. (11) 
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FIG. 6: Linear growth rate A vs wavenumber K from Q 
for different values of fi close to fip. We take ii = 0.1 and 
D = 10"^ so that fip = 0.185 and Km = 47.79. 
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FIG. 7: Spherically averaged continuum structure function 
against K for different values of the control parameter fi. The 
other parameter values are D = 10"^ and R = 0.1. 



Jo and J2 are the zero and second order Bessel functions, 
respectively. This equation can be solved numerically to 
obtain 



so that 

fiP w 185.192^. (13) 

The behavior of the deterministic Eq. {Tj) is now clear: 
for yU < the only stable solution is (p — 0. In the 
interval < fi < iip one has the homogeneous density 
= 0s, and for /i > /ip spatial patterns emerge. This 
last transition can also be crossed by decreasing D/R^ 
at fixed yit > 0. 

This scenario is consistent with the results for the par- 
ticle ND model shown in Sec.|nJ In particular, note that 
Eq. (|12|l indicates that the pattern periodicity is deter- 
mined by i?, and is independent of other parameters of 
the system such as fi, D and Ns- This is in agreement 
with the results for the structure factors shown in Fig. O 
It is also observed there that the numerical value of the 
dominant wavenumber in Fig. is close to the predicted 
value given by Eq. ((T^ . Since Eq. Q has no particu- 
lar symmetries, we expect on generic grounds jig that 
hexagonal patterns would appear close to the instabil- 
ity. Since they usually bifurcate subcritically we expect 
some range of bistability for /i < /ip, that may be in- 
fluenced by noise. In consequence we do not expect the 
transition line (|13|l to be fully accurate. Nevertheless it 
correctly explains the distinct behavior between the data 
shown in Fig. |3t (essentially all of them predicted to be 
in the homogeneous phase, as confirmed by Figs. 01 and 
Isj and Fig.|2|3 (for smaller D, so that all data points are 
in the periodic clustered phase) . The curves in Fig. 
collapse together and approach the deterministic predic- 
tion for the homogeneous solution 0^ (a straight line of 



slope 1 in that scaled plot) sufficiently far from the ab- 
sorbing transition point. Such collapse docs not occur 
in Fig. since they do not correspond to homogeneous 
states. More important are the fluctuation corrections 
to our deterministic results around the absorbing tran- 
sition: the transition point is quite far from the deter- 
ministic value /ic = and the critical behavior is very 
different from the simple linear vanishing of the number 
of particles predicted dctcrministically. 

In Fig. [3 we plot the spherically averaged structure 
function, Sc{K), against the wavenumer, K, of the den- 
sity field (f) obtained numerically, after a long-time, from 
numerical solution of Eq. Sc{K) is the modulus of 
the spatial Fourier transform of (j){:x.,t), averaged spher- 
ically and in time. Note that, since is a continuous 
field, Sc{K) is related but not identical to the structure 
factor S{K) of the particle system, Eq. 10. Neverthe- 
less, maxima of Sc{K) also identify dominant periodic- 
ities. In Fig. [3 we have taken R = 0.1, D — 10~^, so 
that /ip ~ 0.185. One can see how for /t > /tp the 
structure function develops a peak that grows with /i, 
indicating the development of a spatial pattern with a 
typical distance between clusters. The peak is located 
at the wavenumber closest to H12|) compatible with the 
discretisation imposed by the periodic boundary condi- 
tions. Fig. |S1 shows a steady pattern of density which is 
analogous to the one shown for the discrete model in the 
bottom-right panel of Fig. 0] These observations con- 
firm for the full nonlinear model Q the behavior iden- 
tified from the linear stability analysis of the homoge- 
neous solutions. It is also worth to mention that the 
non-sphcrically-averaged version of the structure func- 
tion displays hexagonal order for /t > /ip. It degenerates 
into square symmetry as we increase the value of /i, in- 
dicating the appearance of additional bifurcations. 




FIG. 8: Steady spatial pattern from the deterministic equa- 
tion 0. II ^ 0.70, R = 0.1, D = 1Q~^ and Ns = 50. Note 
the strong similarity with the pattern in the bottom-right plot 
in Fig.El 



V. INFLUENCE OF FLUID FLOW 

In addition to the pattern-forming instability, a cru- 
cial difference between the BB and the neighborhood- 
dependent model is their response to an external flow. 
Since the birth and death rates of the BB model are fixed 
constants, global quantities such as the total number of 
particles are independent of any particle motion, being 
it diffusive or hydrodynamic. This seems to be rather 
unrealistic for applications such as modelling plankton 
populations 2], always driven by external flows. On the 
contrary, with the neighborhood dependence of the rates 
in the tribal model we overpass this inconvenience and 
the model becomes dependent on the environmental con- 
ditions. 

As a simple illustration of the impact of a velocity field, 
we consider the flow given by the Harper map, which is 
a symplectic map in two dimensions and, therefore, it 
resembles an incompressible flow. At each time step the 
fourth item in the algorithm described in Sect. Ullconsists 
in moving the particles in the following way: If we denote 
by {xi{t),yi{t)) the coordinates of the particle i at time 
t, after one iteration of the map they become 



Xi{t') = Xi{t) + Acos{yi{t)), 

y^{t') = y^{t)+AcOs{x^{t')). 



(14) 
(15) 



where t' = t + t. A gives the strength of the flow, 
and depending of its value particles can follow regular 
or stochastic trajectories. Here we are just interested in 
highlighting the behavior when the flow changes, that is 
when A changes. The time evolution of the total number 
of particles is given in Fig. |51 It is seen that the asymp- 
totic value depends on the flow strength A. The reason 
can be clearly seen from Fig. EH where snapshots of the 
particle distributions are presented. It is seen that the 
periodic array of clusters in the absence of flow becomes 
more filamental-like as the flow strength increases. The 
shape of the filamental structures reflects the known un- 
stable and stable foliation of phase space for the Harper 
map. As for the BB model inhomogeneity persists for 
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FIG. 9: N{t), normalized with fiNs/ivR^, vs time for different 
values of the external flow strength, A. From top to bottom: 
A — 0, A = 0.01, A = 0.05, and, fluctuating around the value 
1 (white line), A = 1 and A = 3. The other parameters: 
/i = 0.8, Ns = 50, and R = 0.1. 



rather strong flow, but finally the distribution becomes 
homogeneized. At this point, the particle distribution 
should be very close to Poissonian, with density given by 
the homogeneous solution of 0. This is indeed what 
is observed in the figure (for large values of A the total 
number of particles, N{t), normalized with the survival 
steady density value, fluctuates around 1). For smaller 
flow strength the spatial structure in the neighborhood 
of each particle becomes relevant, and the number of par- 
ticles approaches the value in the absence of flow, given 
in Fig. 13 For completeness. Fig. 1111 shows, in terms of 
the structure factor, the disappearance of structure as 
the flow strength increases. 



VI. SUMMARY 

In this Paper we have introduced an interacting parti- 
cles model which considers birth and death rates for any 
individual depending on the number of particles that are 
within a distance smaller than a given range of interac- 
tion. The model can be considered as a simple non-local 
interacting extension of the BB model, but showing, how- 
ever, a very different behavior. A striking feature is the 
appearence of a quasistationary (fluctuations are present) 
periodic pattern of particle clusters. 

In order to deepen in the understanding of the pat- 
tern forming instabilities and structures in the discrete 
model, we have derived the continuum equation describ- 
ing the particle distribution. Two main features charac- 
terize this equation: a) the presence in the deterministic 
part of an integral term taking into account the non-local 
interaction among the particles, and b) the very complex 
structure of the noise term, which also reflects the non- 
locality of the interaction. To understand the pattern 
forming process we have studied the deterministic part 
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FIG. 10: Long-time spatial structure for the ND model with 
an external flow. From top to bottom and left to right, A = 0, 
A — 0.1, A = 0.5, and yl = 1. The other parameters: fi = 0.9, 
D = 5 X 10"^ Ns = 50, and R = 0.1. 
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FIG. 11: Structure factor for the patterns shown in Fig. 1101 



of the equation. At this level, the interpretation of the 
continuum field as the density field of particles is valid. 
Then, periodic clustering appears as a finite wavelength 
instability in the density equation. 

In addition, we have also considered the response of 
the model to an external driving flow. Global properties 
of the system like the total number of particles depend 
on the flow, via the influence of it in the spatial structure. 

Future extensions of this work will include a further 
study of the instabilities, an explicit consideration of the 
noise term in the continuum description, and the con- 
sideration of different (chemical or biological) species of 
particles interacting through a finite interaction range. 
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APPENDIX: FIELD THEORETICAL APPROACH 



In this Appendix we derive an approximate continuum 
field equation describing the dynamics of our model, by 
using Fock space techniques. We present more results 
than needed for the analysis of Sect ion Hvl but we include 
them here by its intrinsic interest and for future reference. 
In addition to the usual particle interaction and diffusion 
processes, we consider in this Appendix also the process 
of advection, which is not usually addressed within the 
Fock formalism. 

The standard formulation of Fock space methods 
^3 requires the microscopic model to be defined on 
a lattice. Thus we first modify our off-lattice model to 
fit into a grid, apply the standard procedure, and then 
let the grid spacing to become negligibly small, so that 
we recover the continuum ofF-lattice limit. This is some- 
how different from a macroscopic limit, since only the 
grid size, but not the basic lengths in our model, such as 
the interaction range R, will be made to vanish. From 
the lattice approximation we can write the stochastic dy- 
namics in terms of a Hamiltonian operator. A path inte- 
gral representation can then be obtained and a Langevin 
equation can be derived from it. The off-lattice limit is 
then conveniently taken. The whole process can not be 
performed exactly and we need to introduce two approxi- 
mations, one valid for not too large particle densities and 
another restricting to Gaussian fluctuations. 

We divide space in cells forming a lattice of A sites 
and describe the state of the system by the number of 
particles inside each lattice cell {Ni}.i=i^,,,^A- The lattice 
model will be equivalent to the off-lattice one in the limit 
of small lattice size, so that only zero or one particles 
will be present at each site. In addition it is convenient 
to work in continuous time, so that a continuous time 
Markovian Master Equation describes the dynamics. At 
long times (3> t) there should be no difference between 
the discrete and the continuous time approaches. 

Since statistically independent processes are repre- 
sented by additive terms in the Master equation, we can 
treat separately diffusion, advection, and particle trans- 
formations, and then sum up the corresponding contri- 
butions. We begin with the last process: 
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A. Birth and death 

At each site we have a birth rate per particle 



A,: 



max 0, 



and a death rate per particle 



/3i = max 0, /3o 



(16) 



(17) 



The sum is over the cells that are within a distance R 
from the site i (excluding cell i, so that in the continuum 
limit we recover the original prescriptions Q and ^ ). 
Defining P(A^i, TVs, ^a; = P{{N};t) as the proba- 
bility of having A^i bugs in site 1,...,Na bugs in site A at 
time t, we have the following Master Equation: 



dP{{N};t) 
dt 



N 



= ^(iV, - 1)A,(A^, - 1)P(.... 

i=l 
N 

+ ^(iV, + l)f3,{N, + 1)P(..., 

i=l 

N 

- ^A^,A,(iV,)-P(-,^.,-;0 

i=l 
N 

- ^A^,/3(7V,)P(...,iV„...;t) 



■t) 
■t) 



(18) 



The Fock space representation starts by defining the 
many-particle state vector: 

A 

\^>= P{{N};t)ll'^4f'\o>, (19) 

Ni...Na i=l 



where aj and at are the creation and annihilation opera- 
tors at lattice site i, that is, 

al\Ni..M-NA >= \Ni..M + 1-Na >, 
ai\Ni..M.:NA >= N,\Ni...N, - 1-Na >, 



verifying bosonic commutation rules 

[tti, a]] = Sij, 
a,|0 >= 0. 



(20) 
(21) 
(22) 

(23) 
(24) 



Thus we have: 

d\ip> 



dt 



After simple algebra we have a Schrodinger-like equa- 
tion in Euclidean time: 



dt 

with the Hamiltonian: 



-H\i;>, 



(26) 



a ai 



(27) 

Ai and Pi are the operator versions of A^ and /3i, that is, 
expressions p6|l and p7|l with all the particle numbers 
Nk replaced by the number operators a^ak- They can 
be defined from a power series representation of a conve- 
niently regularized version of the non-smooth expressions 
(PI) and (ini). 

To obtain a path integral representation, a first step 
is to use the commutation relations in H until obtaining 
normal ordering (i.e., creation operators to the left and 
annihilation to the right). We call the resulting expres- 
sion Hno- Then an action depending on the classical 
complex variables V'i*(0 ^'^d "02 (0 computed as: 



(28) 

The 'arrow' notation means that the operators a| and 
should be replaced by the indicated classical variables. 
The action allows to calculate average of lattice quanti- 
ties such as Ni as adequate path integrals over ipi and ip* 
involving the weight . For Hamiltonian H27fl . it turns 
out that normal ordering leads to intricate expressions 
in any regularized version of H16|l - H17ll . We note that for 
low enough densities (or large N^), the maximum condi- 
tion would be rarely needed. In consequence, a sensible 
approximation at low density is: 



A,: « An — 



(29) 



and 



jeR{i) 



a ^ I 
— > a c 

TV,- ^ J 



(30) 



{N},t 



In regions where densities are not small, these expressions 
(^]\f^ „ I'jX^l^JSf- — 1)P( Ni — 1 ) Yl^ (a^)^'|0 >would need corrections. We will comment more on this 

*~ ' later. 

The fact that Ai and Pi do not depend on bosonic op- 
erators on site i (the sum j G R(i) excludes the cell i) 
makes trivial the normal-ordering procedure. The action 
(2^ 



~NMNr)Pi...N,...)UtA<)'''\0> 

-im + i)A(iv, + i)p(7v, + i)nti(«I)'^*io > 

-7V,A(^.)P(...iv,...)nti(«l)^'|o>l 
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S = lo dt [^t^^ + (/?o - Ao) 



.jeR{i) v^J 



(31) 



An extra term (not written) should be added to imple- 
ment the particular initial condition used for P{{N},t). 
One could approach the continuum limit at this step 
but, for clarity, we will perform it after writing down 
the Langevin equation. By introducing Gaussian noises 
at each site, {rikit)}, as Gaussian integrals in averages 
involving Ea. H31|l . one realizes that averages of physical 
quantities, such as < Nk > can be obtained as averages 
over stochastic processes {V'fc(^)} satisfying Ito stochas- 
tic differential equations. Terms in the action linear in 
{V'fe} give rise to deterministic terms in the stochastic 
equation, terms quadratic in these variables determine 
the correlations of the Gaussian noises {??fc(t)}. Terms 
of higher order give rise to non-Gaussian noise statistics. 
As usual we neglect these terms (this is our second ap- 
proximation, by which we neglect the last term in (|31|l l 
and obtain an approximate Ito Langevin equation: 



where the Gaussian noises {rjkit)} verify: 

[S^J (Ao - ^ EjeR(i) V'j) 

-%(i?)^^.^,}. 



(32) 



(33) 



(34) 



9ij (i?) is a function with value 1 if cells i and j are differ- 
ent and each one is within a distance R from the other. 
Noise correlations are multiplicative, nonlocal, and rather 
involved. In addition, they can not be satisfied by real 
stochastic process {r]k{t)} so that the noise terms, and 
the variables {t/jkit)}, are in general complex valued. De- 
spite this, the particle statistics is encoded in the lattice 
stochastic process. For example, the probability for the 
occupation numbers is given by 



P{{N},t) 



n 



Ni 



(35) 



where the average is over the statistics of the processes 
{ipi}. Consequences of (p?5ll are < Ni{t) >=< ifjii^t) > 
and < Niit)"^ >=< i/'i(0^ > + < >• The term 

dependent on the initial condition that was omitted from 



(|31|l can now be taken into account by providing adequate 
initial conditions for {^^(t)}, linked to the initial particle 
statistics by (|^ . 

The ofF-lattice continuum limit is performed by intro- 
ducing the density (f>{x,t) by the change tpi 0(x, i)A'^ 
and taking the lattice spacing going to zero A 0. d 
is the spatial dimension {d = 2 through this Paper). We 
have some freedom in performing this limit. For example, 
if we take the interaction range i? to be a fixed number 
of lattice spacings, then i? — > in the continuum limit. 
In this case non-locality is lost, and by properly scaling 
the interaction parameters Ng and a, we get a continuous 
equation related to Reggeon field theory that has been 
thoroughly studied^JE^- But we want to describe the 
macroscopic behavior of the off-lattice models introduced 
in Sect, m in which the interaction range is finite. Thus, 
we fix i? to a finite value when going to the continuum 
and obtain 



N 



with 



(77(x,i)) =0 
(?7(x,i)r,(x',t')) -2<5(t-i')x 



{<5(x-x')0(x,t) [Ao-^/, 



J0<|x-r|<fl 



dr0(r, t) 



(36) 



(37) 



(38) 



The continuum description of the BB model is recov- 
ered if i? ^ 0. As in the discrete case, correlations in H38|l 
are multiplicative, nonlocal, and lead to complex valued 
processes. 

We see that the deterministic part of H36() is what one 
would guess as a mean-field description for the particle 
density when taking into account the finite range of the 
interactions. Since we expect this to be correct for high 
enough densities and far from transition points, we con- 
clude that the terms neglected in the low density approx- 
imation H29l) - (|30l) are just correcting fluctuation statistics 
in a regime, high densities, in which they are already not 
very relevant. 



B. Diffusion and advection 

Advection is implemented as a deterministic process 
in the off-lattice models of Sects. lUll and 0). When 
going to a lattice description it can be implemented, for 
example, as cellular automata deterministic rules. But 
we prefer to model it as a stochastic birth-death process 
in the lattice, that in the continuum limit recovers the 
deterministic flow, since this leads to a more natural de- 
scription in the Fock space formalism we are using, and 
the connection with the diffusion process is clearer. 

By considering stochastic hopping from cell i to j, con- 
taining Ni and Nj particles, at transition rate p{i j) 
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per particle, one is left with the master equation: 

^P(iV„ Nj) = P{N, + 1, Nj - im + IM* ^ j) 

- P{N,,N,)N,p{i^j). (39) 

By considering the time evolution of the Fock space vec- 
tor IV- P[N.,,Nj]t){a\)'^'{a])'^'\0 >, one arrives to 

the Hamiltonian H — ~p{i — s- i)(a] — By consid- 

ering the statistically independent hopping between an 
arbitrary number of cell pairs, one gets 

-Chopping = X! ~ [^(-^ ^ ~ ^ 

<ij> 

where the sum is over all different pairs of sites among 
which there is hopping. In the case in which all the hop- 
ping rates arc equal wc recover the standard form of the 
diffusion Hamiltonian 8, 21] . 

We considerer the case of a Cartesian lattice, with 
cells labelled by the position x and jumps among nearest 
neighbors. In this notation (|40|) becomes 

^hopping = 2 Ex X]iy=l | ("^I+e^A ~ "^x) ^ 

[p(x + e^A ^ x)ax+e„A - p(x ^ x + e,,A)ax]} .(41) 

{Gu}u=i,...,d are unit vectors along the positive direc- 
tions of the cartesian axes. The sum in Eq. H41|) takes 
into account all pairs of nearest neighbors. We split the 
jump rates into a constant part and the local anisotropics 
(which will give rise to the macroscopic flow): 

p(x-^x + e^A) = K + a+(x) (42) 
p(x + A ^ x) = K + (x + e,. A) (43) 

We now perform the continuum limit by letting the lat- 
tice space A ^ 0. To this end we expand the functions 
a^(x -|- e^A) and the operators ax_|_e a ax+e^A in 
powers of A to get 

-^hopping — 

X^x {-^Vat • ((3+(x) - a-(x))ax + A^AtVat • Vcx 
+ ^ ELi V, [(a+(x) + a-(x)) aj 

+0(A3)} . (44) 



d are vectors of components a^. We define the quan- 
tities D = kA^, v{x) = A (a+(x) - a"(x)), and S{x) = 
A^ (q?+(x) + a^(x)) /2, and impose them to be finite in 
the A ^ limit. We apply the recipe H46|l by intro- 
ducing the classical continuum fields a 0(x)A'^ and 
— > (/'*(x) + 1 (Hamiltonian 1)44(1 is already in normal 
order). In the limit A ^ wc find the following expres- 
sion that should be added to the continuum limit of H31() 
to find the complete action of the model: 



f^^^hopping 

JdtJ dx{-V0*(x) • 'y(x)0(x) + D\7(j)*{x) ■ V0(x) 

+ ELi (x)V, (S'.(x)0(x)) } . (45) 

The Langevin representation can now be obtained. All 
the terms in 1)45(1 are of first order in the variables (p* , so 
that diffusion and advection only enter into the determin- 
istic part of the Langevin equation, not in the noise corre- 
lations. In the simpler case in which S'(x) — (for exam- 
ple this happens if a+ (x) = —a~ (x) or if a+ (x) + a~ (x) 
does not diverge fast enough in the A limit) the 
complete Langevin equation reads: 



5t<^(x, <) = -V • (t/(x, t)0(x, t)) + i^V2<^(x, t) 
+ (Ao - /3o) H^, t) 
-^0(x, t) L,^^j,drcbir, t) + 77(x, t), (46) 



which allows the identification of ^(x, t) with a macro- 
scopic velocity flow (we have allowed for an explicit time 
dependence since this does not alter the above deriva- 
tion). The form in which the velocity field appears in 
((46(1 could be guessed from mean- field arguments. The 
important point of the derivation is that it confirms that 
there are no extra terms in the noise correlation arising 
from advection nor diffusion. The noise correlations are 
again given by l(T7|l - l(l^ . 
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